Lower Critical Dimension of Ising Spin Glasses 
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Exact ground states of two-dimensional Ising spin glasses with Gaussian and bimodal (± J) distri- 
butions of the disorder are calculated using a "matching" algorithm, which allows large system sizes 
of up to N = 480 2 spins to be investigated. We study domain walls induced by two rather different 
types of boundary-condition changes, and, in each case, analyze the system-size dependence of an 
appropriately defined "defect energy", which we denote by AE. For Gaussian disorder, we find a 
power-law behavior AE ~ L , with 8 — —0.266(2) and 8 = —0.282(2) for the two types of boundary 
condition changes. These results are in reasonable agreement with each other, allowing for small 
systematic effects. They also agree well with earlier work on smaller sizes. The negative value 
indicates that two dimensions is below the lower critical dimension d c . For the ±J model, we obtain 
a different result, namely the domain- wall energy saturates at a nonzero value for L — > oo, so 8 = 0, 
indicating that the lower critical dimension for the ±J model exactly d c = 2. 



PACS numbers: 75.50.Lk, 05.70.Jk, 75.40. Mg, 77.80.Bh 



Although most of the effort to understand spin glasses 
PI has been concerned with the behavior in three dimen- 
sions, the situation in two dimensions is also not fully 
understood. For a Gaussian distribution of disorder, it is 
widely accepted that g-Q], T c — and two-dimensions 
is below the lower critical dimension d c . However, for 
bimodal (±J) disorder, inconsistent results have been 
found. Using exact ground-state calculations of small 
systems (N < 30 ) and introducing domain walls in the 
system by changing the boundary conditions, Kawashima 
and Rieger Q also concluded that T c = and probably 
d c > 2, though they could not rule out marginal behav- 
ior, i.e. d c = 2. On the other hand, from Monte-Carlo 
(MC) simulations || on sizes N < 18 2 , it was claimed 
that T c > 0. A calculation of domain-wall energies us- 
ing a cluster MC method J7J for sizes N < 24 2 indicated 
marginal behavior (i.e. d c = 2), and this was interpreted 
as evidence for the finite- T transition claimed in Ref . || . 
Recently, Houdayer || used a more efficient cluster MC 
algorithm, which allows larger systems (N < 100 2 ) to 
be studied down to low temperature T = 0.1. Houdayer 
finds that T c = but with an exponential divergence 
of the correlation length for T — » 0, which is consistent 
with d c = 2, in agreement with earlier work by Saul and 
Kardar ||. 

To try to clarify the situation, we calculate here 
ground-state domain-wall energies AE using an exact 
polynomial-time matching algorithm |L2] for much larger 
system sizes (N < 480 2 ) than in previous work. For the 
Gaussian model our results are consistent with earlier 
calculations and we find that T c = and d c > 2. For the 
±J model, we find marginal behavior expected at the 
lower critical dimension, and so conclude that d c = 2. 

The model consists of N — L 2 Ising spins Si = ±1 on 
a square lattice with the Hamiltonian 



ri — 2_^ JijSiSj , 
(i,3) 



(1) 



(i,j) and the Jij are quenched random variables. Here, 
we consider two kinds of disorder distributions: (i) Gaus- 
sian with zero mean and variance one, and (ii) a bimodal 
distribution J,^ = ±1 with equal probability. 

In greater than two dimensions, or in the presence of a 
magnetic field, the exact calculation of spin-glass ground 
states belongs to the class of NP-hard problems |10|.|llJ . 
This means that only algorithms with exponentially in- 
creasing running time are known. However, for the spe- 
cial case of a planar system without magnetic field, e.g. 
a square lattice with periodic boundary conditions in at 
most one direction, there are efficient polynomial-time 
"matching" algorithms p2|. The basic idea is to rep- 
resent each realization of the disorder by its frustrated 
plaquettes [^3|. Pairs of frustrated plaquettes are con- 
nected by paths in the lattice and the weight of a path 
is defined by the sum of the absolute values of the cou- 
pling constants which are crossed by the path. A ground 
state corresponds the set of paths with minimum total 
weight, such that each frustrated plaquette is connected 
to exactly one other frustrated plaquette. This is called 
a minimum- weight perfect matching. The bonds which 
are crossed by paths connecting the frustrated plaquettes 
are unsatisfied in the ground state, and all other bonds 
are satisfied. 

For the calculation of the minimum-weight perfect 
matching, efficient polynomial-time algorithms are avail- 
able p3Jq| . Recently, an implementation has been pre- 
sentedfl^f, where ground-state energies of large systems 
of size N < 1800 were calculated. Here, an algorithm 
from the LEDA library flTJ has been applied, which lim- 
its the system sizes to N < 480 2 , due to the restricted 
size of the main memory of the computers which we used. 

To study whether an ordered phase is stable at finite 
temperatures, the following procedure is usually applied 
||| |,^8| |20). First a ground state of the system is calcu- 
lated, having energy E . Then the system is perturbed to 
introduce a domain wall and the new ground state energy, 

pertb 



where the sum runs over all pairs of nearest neighbors E ° is evaluated. Typically, the system initially has 




FIG. 1. A domain wall created in an L = 320 system with 
Gaussian interactions. Free boundary conditions are applied 
across the top and bottom edges (dashed lines) and periodic 
boundary conditions applied across the vertical edges (solid 
lines). The domain wall is induced by changing the periodic 
boundary conditions to antiperiodic, which is equivalent to 
changing the sign of the bonds which wrap round the system 
from left to right. 



periodic boundary conditions in both directions, and the 
perturbation involves replacing periodic by anti-periodic 
boundary conditions in one direction. However we can- 
not apply the matching algorithm for boundary condi- 
tions which "wrap around" in both directions, so instead 
we have investigated two slightly different procedures: 

(i) We choose boundary conditions which are free in 
one direction and periodic in the other (x-say). 
The periodic boundary conditions are then replaced 
by anti-periodic, which is equivalent to changing 
sign of the bonds which "wrap around" the sys- 
tem. This is very similar to the usual approach, 
except that the induced domain wall no longer has 
to have the same x coordinate at the top and bot- 
tom of the sample. An example of a domain wall 
formed in this way is shown in Fig. [I]. We denote 
these as "P-AP" boundary condition changes. Note 
that the energy change can have either sign, and 
the distribution of values (obtained by repeating 
the calculation for many samples) is symmetric. 

(ii) We apply free boundary conditions in both direc- 
tions and compute the ground state. We then add 
extra bonds which wrap around the system in the 
x-direction say, and have a sign and strength such 
that they force the spins they connect to have the 



opposite relative orientation to that which they had 
in the original ground state. The new ground state 
is then computed and we remove from the ground 
state energy the contribution from the added strong 
bonds. The change in energy must then be posi- 
tive. This procedure is similar to the "conjugate 
boundary conditions" chosen in Ref. M . Since free 
boundary conditions are still applied along the top 
and bottom edges, the domain walls look very sim- 
ilar to the one shown in Fig. |]. We denote this 
boundary condition change by "F-DW", since the 
initial boundary condition is Free, and it is then 
changed to induce a Domain Wall. 

For a given sample, the domain- wall energy is given by 

SE = £P° rtb - E . (2) 

To study finite-size behavior we consider the "defect en- 
ergy" , AE, defined by 



AE=[\SE\] a 



(3) 



where [• • -] av denotes an average over samples. (Of course 
the absolute value is not necessary for the F-DW bound- 
ary conditions since SE > 0.) It is expected that AE will 
vary as 



AE ~ L\ 



(4) 



where 9 is a stiffness exponent. For 9 < an ordered 
phase is only stable at T — and the correlation length 
£ diverges g for T -> as £ - T 1 / 6 '. For 9 > an 
ordered phase also occurs for T > 0. At the lower critical 
dimension, one has 9 = and expects an exponential 
divergence of the correlation length. 

Domain-wall energies obtained in the way described 
above were obtained for systems with Gaussian and ± J 
disorder, for both F-DW and P-AP boundary condition 
changes. Sizes up to N — 480 2 were considered, and, 
for each size, the result was average over up to 40000 
samples. Table B shows the number of samples for each 
system size and type of systenr. 



type 


L < 160 


L 


= 240 


L 


= 320 


i = 480 


±J P-AP 


40000 




45000 




24000 


5000 


±J F-DW 


10000 




5500 




4500 


- 


Gaussian P-AP 


10000 




2800 




3700 


- 


Gaussian F-DW 


10000 




700 




- 


- 



TABLE I. The number of samples for each type of disorder 
and system size. 

In Fig. g the defect energies are shown as a function of 
system size. We first discuss the results for the Gaussian 
distribution. The domain-wall energy clearly decreases 
with a power law AE ~ L e with 



9 = -0.282 ±0.002 
9 = -0.266 ±0.002 



(P - AP) 

(F-DW). 



(5) 
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FIG. 2. The defect energy AE as a function of system 
size for the different types of disorder and boundary condi- 
tion changes. The lines are fits (for L > 8) to the form 
AE = AEoc + bL . The following fit parameters were ob- 
tained: AE = 0.96, b = 0.99, 9 = -0.65 (±J P-AP); 
AE = 1.66, b = 2.93, 9 = -0.99 (±J F-DW); b = 1.58, 
9 = -0.282 (AE = 0, Gaussian P-AP); b = 2.78, 9 = -0.266 
(AE = 0, Gaussian F-DW). 



For the goodness of the fits |2l| we obtain Q = 0.48 and 
Q = 0.62 for the P-AP and F-DW boundary conditions 
respectively. In addition, we tried fits of the form AE — 
A^oo + bL e , obtaining AE X = 0.04 ± 0.02 (Q = 0.48) 
for P-AP boundary conditions and AE^ = -0.03 ±0.07 
(Q = 0.54) for F-DW boundary conditions. Hence, we 
cannot rule out that the domain-wall energy converges to 
a small nonzero value for L — > oo but this seems unlikely. 

The result of the exponent for the P-AP boundary con- 
ditions is compatible with the value of = —0.281(2) 
found H for full periodic boundary conditions and sys- 
tem sizes L < 30. The value of the F-DW case disagrees 
with the value for the P-AP case by more than the er- 
ror bars, but the difference is quite small. We expect 
that this is due to systematic corrections to scaling and 
that asymptotically the same result would be obtained, 
i.e. the same exponent would be obtained for any rea- 
sonable definition of the defect energy. Our result for 
the F-DW boundary condition differs significantly from 
the value 6 = —0.20 obtained J7J]. The reason for the 
discrepancy may be that Ref. jj studied smaller sizes, 
L < 24, and used a Monte-Carlo method to determine 
the ground states instead of an exact algorithm. Nev- 
ertheless, to our knowledge, all results obtained for the 
Gaussian model are compatible with T c = 0. 

Next we turn to the results for the ± J model. The dif- 



ference between these results and those for the Gaussian 
model, both of which are shown in Fig. |2| is quite strik- 
ing. Whereas the data for the Gaussian model clearly 
tends to zero with a power law for L — ► oo, the data for 
the ±J distribution appears to tend to a non-zero value. 
Fits of the same form as for the Gaussian model give 
AEoc = 0.96 ± 0.01 (Q = 0.23) for the P-AP boundary 
conditions and AE^ = 1.66 ±0.01 (Q = 0.50) for the F- 
DW boundary conditions. Fits in which AE^o was fixed 
to zero were also tried but gave ridiculously low probabil- 
ities: Q = 2 x 10~ 27 and Q = 2 x 10~ 17 respectively. For 
the case with free boundary conditions in the x-direction, 
a saturation has been observed before [R for small sys- 
tem sizes L < 24, using a Monte-Carlo method to find 
ground states, rather than an exact algorithm. 
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FIG. 3. Distribution P(AE) of domain-wall energies for 
the ±J distribution with F-DW boundary conditions. The 
points for L — 320 and L — 10 are slightly displaced on the 
x-axis for visibility. 

To show the effect of saturation into more detail, Fig. || 
displays the distribution of the domain- wall energies, 
with a logarithmic vertical scale, for the ±J distribution 
with F-DW boundary conditions. For this distribution, 
only even integer values of AE are found. Except for the 
largest values of AE, where the weights are very small, 
the weights of the discrete bins are almost the same for 
L — 80 and 320 indicating that the whole distribution 
converges for large L. The probability P(0) that the 
domain-wall energy is exactly zero is displayed in Fig. H. 
P(0) seems to converge towards a finite value less than 
unity for both boundary conditions. 

From this data we deduce that 9 = for the ± J distri- 
bution, so the lower critical dimension is d c = 2 exactly. 
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FIG. 4. The probability P(0) that the domain-wall energy 
is exactly zero for the ± J model as a function of system size. 
The lines are guides to the eye. 



The most likely scenario at finite temperature is that the 
correlation length diverges exponentially as T — > 0, per- 
haps like exp(C/T) or exp(C/T 2 ). Detailed Monte Carlo 
simulations g and the work of Saul and Kardar || indi- 
cate that the former is the case. The other possibility is 
a finitc-T Kostcrlitz-Thouless type of transition without 
true long range order below T c , as claimed by Ref. ||. 
However, the finite value for v obtained in Ref. M seems 
incompatible with this picture. The reason for this dis- 
crepancy is probably that the system sizes were quite 
small, N < 18 2 . 

To conclude, we have studied the defect energy of the 
two-dimensional Ising spin glass using much larger sys- 
tem sizes than before, up to N = 480 2 . For Gaussian 
disorder we find that a spin-glass phase is not stable 
at T > 0, and that d c > 2, in agreement with earlier 
findings for much smaller systems. For the ±J model, 
the absolute value of the domain-wall energy saturates 
at a nonzero value for L — > oo, indicating that the lower 
critical dimension is exactly d c = 2, in agreement with 
Refs. B and |9j. It is quite striking and perhaps surpris- 
ing that the ±J model, which has a large ground state 
degeneracy, actually has more order at very low-T than 
the Gaussian distribution, which has a unique ground 
state. 
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